* Gary LaMotte's Mathematica solutions to exercises from the 2012 Gravitational Wave Astronomy Summer School
* Day 1 part 3

* WELCH'S METHOD NOISE + SINES*

*RESULTS OF THIS EFFORT SHOW THAT BY SEGMENTATION AND OVERLAP WELCH'S METHOD DOES NOT CHANGE THE MAGNITUDE OF THE NOISE AS MUCH AS EXPECTED, BUT RATHER FLATTENS THE NOISE FLUCTUATIONS WHILE LEAVING THE PEAK HEIGHT OF THE SIGNAL LITTLE CHANGED.  THE PEAK IS BROADENED SUBSTANTIALLY.*

samprate = 16384;

nsec=1.0;

ndata=nsec*samprate;

time = Table[(r-1.0)/ndata  , {r,ndata}];

sines = 5*Sin[2*Pi*74.7*time] + 2*Sin[2*Pi*820*time];

noise = Table[RandomVariate[NormalDistribution[0,1.4],ndata]];

ts = noise + sines;

st = Transpose[{time, ts}];

ListLinePlot[st, PlotLabel -> Timeseries, AxesLabel ->  {Time (Secs), Amplitude}]

ListLinePlot[st, PlotLabel -> Timeseries, AxesLabel ->  {Time (Secs), Amplitude}, PlotRange-> {{0,0.05}, {-11,11}}]

ftsines = Fourier[ts];

ListLinePlot[(Abs[ftsines])^2,PlotRange -> {{0,16384},{0, 80000}}]

ListLogPlot[(Abs[(ftsines)])^2,PlotRange->{{0,16384},{5,300000}},Joined->True, PlotLabel-> Periodogram, AxesLabel-> {Frequency, Amp^2}]
*Very tedious and long to develop this input.  Seems it should have been much more simple.*

ListLogPlot[(Abs[(ftsines)])^2 ,PlotRange -> {{0,1000},{2, 300000}},Joined->True, PlotLabel-> Periodogram, AxesLabel-> {Frequency, Amp^2}]

*HAMMING FILTER*

filt = (0.54 - 0.46(Cos[2*Pi*time]))*ts;

graph = Transpose[{time, filt}];

ListLinePlot[graph, PlotLabel -> Timeseries, AxesLabel ->  {Time (Secs), Amplitude}]

ftfilt = Fourier[filt];

nyqpart = ftfilt[[1;;8192]];

powerspectrum = (Abs[nyqpart])^2;

ListLogPlot[powerspectrum,PlotRange -> {{0,8192},{.01, 100000}}, Joined->True, PlotLabel-> NonSegmented Hamming]

ListLogPlot[powerspectrum,PlotRange -> {{0,1024},{.01, 100000}}, Joined->True, PlotLabel-> NonSegmented Hamming]

ListLogPlot[powerspectrum,PlotRange -> {{0,130},{.01, 100000}}, Joined->True, PlotLabel-> NonSegmented Hamming]

*BELOW: POWER SPECTRA WITH WELCH METHOD
PURPOSE:  DECREASE NOISE WITH UNFORTUNATE SPECTRAL LEAKAGE OF PEAK*

*Set up 8PARTITIONS and new time length*

partitions = 8;

segments = Partition[ts,     Length[ts]/partitions,      Length[ts]/(2*partitions)];

nsec2 = nsec / partitions

ndata2 =  (samprate * nsec) / partitions

newtime = Table[(r-1.0)*(nsec2/ndata2)  , {r,ndata2}];

*DO HAMMING FILTERING OF EACH SEGMENT BELOW*

filtsegs = Transpose[(0.54 - 0.46(Cos[2*Pi* newtime/nsec2]))*Transpose[segments]];

*Do Hann filtering instead if desired: filtsegs = Transpose[Transpose[segments] * (1-Cos[2*Pi*newtime/nsec2])];*

f = filtsegs[[5]];

graph2 = Transpose[{newtime, f}];

ListLinePlot[graph2, PlotLabel -> Timeseries, AxesLabel ->  {Time (Secs), Amplitude}]

ftfiltsegs = (Abs[Fourier[filtsegs]])^2;

nyqtrans = Transpose[ftfiltsegs][[1;;(samprate/(partitions*2))]];

nyqpart = Transpose[nyqtrans];

added= Total[nyqpart];

averaged = added/((partitions*2) - 1);

index = Table[(r-1.0)*partitions  , {r,Length[averaged]}];

normalized = Transpose[{index, averaged}];

ListLogPlot[normalized,PlotRange -> {{0,8192},{.01, 100000}}, Joined->True, PlotLabel-> SegmentedHamming PSD Welch]

ListLogPlot[normalized,PlotRange -> {{0,1024},{.01, 100000}}, Joined->True, PlotLabel-> SegmentedHamming PSD Welch]

ListLogPlot[normalized,PlotRange -> {{0,130},{.01, 100000}}, Joined->True, PlotLabel-> SegmentedHamming PSD Welch]

*TOTAL OF 8 PARTITIONS IN THE TIME SERIES IN THE ABOVE ANALYSIS*

*COMPARE TO SEGMENTED HANN  WELCH PSD BELOW
NOTE:  PEAKS BELOW HAVE  MORE SPECTRAL LEAKAGE AND MORE NOISE VS SEGMENTED HAMMING WELCH PSD*